clc;
clear;

format short

pho1 = 3;
pho2 = 1;
L = 1;
d = 2;
h = 0.075;

[x, y] = meshgrid(-L:h:L, -L:h:L);

eta = @(x) -0.1*d*cos(2*pi*x./d);

pho = @(x,y) (pho1 + pho2)/2 + (pho1 - pho2)/2*tanh((y-eta(x))/0.01*d);

surf(x, y, pho(x,y));
